%%%%%%%%%%%%%%%%%%
%% SETUP THE MATLAB PATHS
% make sure that fieldtrip and spm are not in your matlab path

% SET THE BELOW LINE TO THE OSL DIRECTORY:
%tilde='/home/mwoolrich';
tilde='/Users/woolrich';

osldir=[tilde '/homedir/matlab/osl1.1'];
addpath(osldir);
osl_startup(osldir);

%%%%%%%%%%%%%%%%%%
%% INITIALISE GLOBAL SETTINGS FOR THIS ANALYSIS
% This specifies where our data is stored, what the data filenames are

datadir=['/home/hluckhoo/OSL_development_scripts/ICA_workshop/faces_subject1_data/fifs']; % directory where the data is

workingdir=[datadir '/spmfiles']; % this is the directory the SPM files will be stored in

cmd = ['mkdir ' workingdir]; unix(cmd); % make dir to put the results in

% Set up the list of subjects and their structural scans for the analysi 
% Currently there is only 1 subject.
clear fif_files spm_files structural_files

% list of fif files (we only have one here - note that it has already been
% Maxfiltered and downsampled to 250Hz)
%fif_files{1}=[datadir '/fifs/sub1_face']; 
fif_files{1}=[datadir '/sub2_rest']; 
% set up a list of SPM MEEG object file names (we only have one here)
spm_files{1}=[workingdir '/spm8_meg1'];
spm_files_sss{1}=[workingdir '/spm8_meg1_sss'];
spm_files_sss_RBC{1}=[workingdir '/spm8_meg1_sss_RBC']; % RBC = Rejected Bad Channels manually
% structural files in the same order as spm_files or fif_files:
structural_files{1}=[datadir '/structurals/struct1.nii'];

cleanup_files=0; % flag to indicate that you want to clean up files that are no longer needed

%%%%%%%%%%%%%%%%%%

%% 1.) Run Maxfilter without rejecting bad channels and Convert Output to SPM.
% This is the wrong way to approach your data. If you use this script as a
% template then you don't want to include this cell. 

for subnum=1:length(fif_files)
Smf=[];
Smf.fif=[fif_files{subnum}];
Smf.logfile=1;
Smf.downsample_factor=4;
fif_sss=osl_call_maxfilter(Smf);
end

S2=[];
for subnum=1:length(fif_files), % iterates over subjects
    S2.fif_file=[fif_files{subnum} '_sss.fif'];
    S2.spm_file=spm_files_sss{subnum};
    
    % The conversion to SPM will show a histogram of the event codes
    % and correspond to those listed below in the epoching section
    [D spm_files_sss{subnum}] = osl_convert_script(S2);
end
close all;

% Delete the bad fif file as we don't need it and later MaxFilter Calls
% will want to write using the same file name.
unix(['rm ' fif_sss '.fif'])

%%%%%%%%%%%%%%%%%%
%% 2.a) Create 1st round of un-maxfiltered data
% In order to look for scanner artefacts we need to generate an
% un-Maxfiltered data set first. 

for subnum=1:length(fif_files)
Smf=[];
Smf.fif=[fif_files{subnum}];
Smf.logfile=1;
Smf.downsample_factor=4;
Smf.nosss=1;
fif_sss=osl_call_maxfilter(Smf);
end

%% 2.b) CONVERT FROM FIF TO AN SPM MEEG OBJECT:

if(length(fif_files)>0),
    S2=[];
    for subnum=1:length(fif_files), % iterates over subjects
        S2.fif_file=[fif_files{subnum} '_nosss.fif'];
        S2.spm_file=spm_files{subnum};       
        [D spm_files{subnum}] = osl_convert_script(S2);
    end;
end;

%% 2.c) Use sensor-space ICA to find the scanner artefacts
for subnum=1:length(fif_files),
    S=[];
    S.fname=spm_files{subnum};
    S.do_plots=1;
    S.do_mains=0;
    S.do_ecg=0;
    S.do_blinks=0;    
    [spm_files_new{subnum} scanner_artefact{subnum}]=osl_africa(S);
end

for subnum=1:length(fif_files),
    D=spm_eeg_load(spm_files{subnum});
    figure;
    for p=1:length(scanner_artefact{subnum})
        subplot(ceil(length(scanner_artefact{subnum})/2),2,p);
        plot(D.time,D(scanner_artefact{subnum}(p),:));
        title(['Scanner artefact ' num2str(p) ' Channel ' num2str(scanner_artefact{subnum}(p))]);
        xlabel('Time (s)');
    end
    selchan=input('Please flag bad channels for removal (e.g. [1 3 5]):     ');
    if ~isempty(selchan)
        chans_marked_as_bad=scanner_artefact{subnum}(selchan)';
        D=badchannels(D, D.meegchannels, ismember(D.meegchannels,[D.badchannels chans_marked_as_bad]));
        %disp(['Channels: ' D.chanlabels(find(ismember(D.meegchannels,[D.badchannels chans_marked_as_bad]))) ' are marked as bad'])
        D.save
    end
end

%% 2.d) Manual Inspection of Channels
for subnum=1:length(fif_files)
    D=spm_eeg_load(spm_files{subnum});
    figure; imagesc(D(1:306,:,:));
    osl_badchans(spm_files{subnum})
    clear D
    D=spm_eeg_load(spm_files{subnum}); D.badchannels
    close all;
end

%%%%%%%%%%%%%%%%%%%
%% 3.) Run Maxfilter, rejecting bad channels and Convert Output to SPM.

for subnum=1:length(fif_files)
Smf=[];
Smf.fif=[fif_files{subnum}];
Smf.logfile=1;
Smf.downsample_factor=4;
Smf.spmfile=spm_files{subnum};
fif_sss=osl_call_maxfilter(Smf);
end

S2=[];
for subnum=1:length(fif_files), % iterates over subjects
    S2.fif_file=[fif_files{subnum} '_sss.fif'];
    S2.spm_file=spm_files_sss_RBC{subnum};
    [D spm_files_sss_RBC{subnum}] = osl_convert_script(S2);
end
close all;

%%%%%%%%%%%%%%%%%%%
%% 4.) Compare the data from stages 1, 2 & 3

Dnosss=spm_eeg_load(spm_files{1});
Dsss = spm_eeg_load(spm_files_sss{1});
Drbc = spm_eeg_load(spm_files_sss_RBC{1});


figure; 

subplot(2,3,1);
mags=find(strcmp(Dnosss.chantype,'MEGMAG'));
[~,order]=sort(std(Dnosss(mags,:),[],2),'descend');
plot(downsample(Dnosss.time,5), downsample(Dnosss(mags(order(1:5)),:)',5)');
xlabel('Time (s)'); title('5 Highest Variance Magnetometers - No SSS')

subplot(2,3,4);
plans=find(strcmp(Dnosss.chantype,'MEGPLANAR'));
[~,order]=sort(std(Dnosss(plans,:),[],2),'descend');
plot(downsample(Dnosss.time,5), downsample(Dnosss(plans(order(1:5)),:)',5)');
xlabel('Time (s)'); title('5 Highest Variance Planar Gradiometers - No SSS')

subplot(2,3,2);
mags=find(strcmp(Dsss.chantype,'MEGMAG'));
[~,order]=sort(std(Dsss(mags,:),[],2),'descend');
plot(downsample(Dsss.time,5), downsample(Dsss(mags(order(1:5)),:)',5)');
xlabel('Time (s)'); title('5 Highest Variance Magnetometers - SSS all channels')

subplot(2,3,5);
plans=find(strcmp(Dsss.chantype,'MEGPLANAR'));
[~,order]=sort(std(Dsss(plans,:),[],2),'descend');
plot(downsample(Dsss.time,5), downsample(Dsss(plans(order(1:5)),:)',5)');
xlabel('Time (s)'); title('5 Highest Variance Planar Gradiometers - SSS all channels')

subplot(2,3,3);
mags=find(strcmp(Drbc.chantype,'MEGMAG'));
[~,order]=sort(std(Drbc(mags,:),[],2),'descend');
plot(downsample(Drbc.time,5), downsample(Drbc(mags(order(1:5)),:)',5)');
xlabel('Time (s)'); title('5 Highest Variance Magnetometers - SSS good channels only')

subplot(2,3,6);
plans=find(strcmp(Drbc.chantype,'MEGPLANAR'));
[~,order]=sort(std(Drbc(plans,:),[],2),'descend');
plot(downsample(Drbc.time,5), downsample(Drbc(plans(order(1:5)),:)',5)');
xlabel('Time (s)'); title('5 Highest Variance Planar Gradiometers - SSS good channels only')

